home *** CD-ROM | disk | FTP | other *** search
/ Monster Media 1996 #15 / Monster Media Number 15 (Monster Media)(July 1996).ISO / prog_c / cuj0696.zip / DWYER.ZIP / SPECTRAL.TST / EUCLSTEP.CX < prev    next >
Text File  |  1996-01-29  |  2KB  |  53 lines
  1. /* ============ */
  2. /* euclstep.cx    */
  3. /* ============ */
  4. #include <xtendefs.h>
  5. #define    P(S)
  6. /* ==================================================================== */
  7. /* EuclStep - Performs Step S2 (Euclidean Step) of the Spectral Test    */
  8. /* ==================================================================== */
  9. void
  10. EuclStep(USHORT *XH, USHORT *XHP, USHORT *XP, USHORT *XPP,
  11.     USHORT *XS, USHORT *XU, USHORT *XV)
  12. {
  13.     for (;;)
  14.     {
  15.     USHORT    T1[NE], T2[NE];
  16.  
  17.     XDIV(XHP, XH, T1);        /* q = h'/ h    */
  18.     XFLOOR(T1, T1);            /* q = floor(q) */
  19.  
  20.     XMULT(XH, T1, T2);        /* h * q    */
  21.     XSUB(XHP, T2, XU);        /* u = h' - qh    */
  22.  
  23.     XMULT(XP, T1, T2);        /* p * q    */
  24.     XSUB(XPP, T2, XV);        /* v = p' - qp    */
  25.  
  26.     XSQR(XU, T1);            /* u * u    */
  27.     XSQR(XV, T2);            /* v * v    */
  28.     XADD(T1, T2, T1);        /* u^2 + v^2    */
  29.  
  30.     if (XGTE(T1, XS))        /* if (s <= u^2 + v^2)    */
  31.     {
  32.         break;            /* Algorithm complete    */
  33.     }
  34.  
  35.     XCOPY(T1, XS);            /* s <- u^2+v^2 */
  36.     XCOPY(XH, XHP);            /* h' <- h    */
  37.     XCOPY(XU, XH);            /* h  <- u    */
  38.     XCOPY(XP, XPP);            /* p' <- p    */
  39.     XCOPY(XV, XP);            /* p  <- v    */
  40.     P(XPRINTF("Next EuclStep: XH = ", XH, NDEC));
  41.     P(XPRINTF("              XHP = ", XHP, NDEC));
  42.     P(XPRINTF("               XP = ", XP, NDEC));
  43.     P(XPRINTF("              XPP = ", XPP, NDEC));
  44.     P(XPRINTF("               XS = ", XS, NDEC));
  45.     P(XPRINTF("               XU = ", XU, NDEC));
  46.     P(XPRINTF("               XV = ", XV, NDEC));
  47.     }
  48.     P(exit(1));
  49.     P(printf("\nExit Eucl: XH = %.Lf, XHP = %.Lf, XP = %.Lf, XPP = %.Lf\n",
  50.         *XH, *XHP, *XP, *XPP));
  51.     P(printf("\t      S = %.Lf,  U = %.Lf,  V = %.Lf\n", *S, *U, *V));
  52. }
  53.